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Abstract 

In this paper we discuss a simple method of testing for the presence of energy-dependent dis- 
persion in high energy data-sets. It uses the minimisation of the Kolmogorov distance between 
the cumulative distribution of two probability functions as the statistical metric to estimate the 
magnitude of any spectral dispersion within transient features in a light-curve and we also show 
that it performs well in the presence of modest energy resolutions (~ 20%) typical of gamma-ray 
observations. After presenting the method in detail we apply it to a parameterised simulated 
lightcurve based on the extreme VHE gamma-ray flare of PKS 2155-304 observed with H.E.S.S. 
in 2006, in order to illustrate its potential through the concrete example of setting constraints 
on quantum-gravity induced Lorentz invariance violation (LIV) effects. We obtain comparable 
limits to those of the most advanced techniques used in LIV searches applied to similar datasets, 
but the present method has the advantage of being particularly straightforward to use. Whilst 
the development of the method was motivated by LIV searches, it is also applicable to other as- 
trophysical situations where energy-dependent dispersion is expected, such as spectral lags from 
the acceleration and cooling of particles in relativistic outflows. 
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1. Introduction 

The timing properties of a sequence of events can be very revealing as to the physical nature 
either of the emitting source or of the medium they propagate through, especially when taken in 
conjunction with information about their energy. Timing analysis algorithms with the capacity 
of resolving energy-dependent properties can then be an important tool for probing the physical 
mechanisms leading to flux variability, such as particle acceleration and cooling. Methods are 
traditionally based on cross-correlation of the binned time-series (e.g. lUt]), and sometimes rely 
on a particular parameterisation of the light-curve, for ex ample by modeling the data according to 
a pre-determined choice for the light-curve profile (e.g. [2]). In the case of gamma-ray sources, 
where high-energy processes are responsible for extreme and short-lived variability events, and 
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for which the observational data are often limited by low photon statistics, unbinned methods are 
the natural and preferential choice of approach to the problem of temporal analysis of time- and 
energy-stamped photon lists. 

In this paper we present a method to search for energy-dependent dispersion in light-curves 
of transient sources with relatively sparse events, particularly suited for (though of course not 
limited to) ground-based gamma-ray telescopes. A particular motivation for the study of energy 
dependent dispersion in the very-high energy (VHE) regime is the prospect of testing for possible 
signatures of Lorentz invariance violation (LIV), foreseen by a number of theories of quantum 
gravity (QG) The analysis method is described in section |2] and its suitability for the par- 
ticularly challenging task of searching for LIV effects are illustrated in section |3] and applied to 
AGN lightcurves in section |4] 

1.1. Lorentz invariance violation 

The unification of the theories of quantum mechanics, governing the smallest of scales, and 
that of gravity, governing the largest of scales, is one of the most serious challenges in modem 
physics. Because of the extremely high energies at which QG effects are expected to manifest 
(around the Planck scale, Eqq ^ By ^ lO'^ GeV) the effects are only likely to become noticeable 
at very high energies and so difficult to be assessed directly in the laboratory. The so-called time- 
of-flight experiments, first proposed in a seminal paper by Amelino-Camelia et al. (1998) |4], is 
one of the most promising ways of carrying out tests for QG signatures. The method is based on 
the search for an energy-dependent speed of light (c) in vacuu from the observation of GeV- 
TeV photons propagating over cosmological distances. The exact form of the energy-dependent 
photon momentum due to QG effects can vary depending on the particular theory adopted, but 
given that its effect is very small it can be treated perturbatively leading to a form (eg, following 
the scheme ofldSiH) 

= £'[1 + ^lEylEgc + UE'ylElc) + ■ ■ ■] (1) 

The consequently small magnitude of its signature at astrophy sic ally accessible energy ranged 
mean that these searches require extremely sensitive measurements. In time-of-flight experi- 
ments the cumulative temporal effects of small variations in c is amplified, eventually manifesting 
as measurable time-delays over the integrated distance travelled by the photons. 

To first order, the magnitude of the time delays expected from QG variations of c are 6t oc 
Ey/EQG ~ 10 s/TeV/Gpc for Planck scale QG. This implies that searches from distant sources are 
preferred (which in turn can lead to them being correspondingly fainter than nearby sources) and 
that the searches should be conducted over narrow features (see section l3T2l for further details). 
For instance, in the case of the active galactic nucleus PKS 2155-304 located at a redshift z ~ 
0.1 16, for which we would expect a delay of 6t ~ 4s per TeV in photon energy, we would need 
flare features on timescales of no more than tens to hundreds of seconds in the VHE light-curve 
to bring the effect to the fore. With event rates of a few Hz during the brightest flares |7] the 
latter property disfavours binning methods on count-rate limited datasets. Sensitivity to small 
spectral dispersions within very limited photon lists is therefore the most desired characteristic 
of a dispersion-search method used for time-of-flight measurements. 



^This is because in QG tlieories the vacuum is expected to have a non-trivial refractive index due to fluctuations of 
the space-time at the quantum level. 

^The most energetic photons recorded from astrophy sical sources have energies of ~ tens of TeV and for Ey ~ I TeV 
the connection to the speed of light due to quantum gravity would be of order 10""c 
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2. Unbinned Methods - Dispersion Cancellation Algorithm 



Unbinned algorithms are well-suited for the identification and analysis of local and aperiodic 
light-curve features, such as bursts or flares in AGN or GRB data. Indeed, the observation 
of GeV photons from GRB 0809 16C and GRB 090510 by the Fermi/LAT coflaboration l^ll 
has recently been able to set limits at and just above the Planck scale for linear-term effects 
using two different unbinned approaches. The first of these methods was to directly compare 
the arrival time of the highest energy photon to the different burst features and, assuming they 
were contemporaneous at the source, determine what magnitude of dispersion would have to be 
experienced by the photons during propagation to the observer to explain the observed lag. There 
are two main drawbacks in such an approach. The first caveat comes from uncertainty in the 
knowledge of the intrinsic structure of the light-curve, due for example to a lack of understanding 
of precursor activity in GRBs (see, e.g. |9]), which can cast doubt as to which particular features 
to associate with the highest energy photons upon assigning the delay. The second drawback 
has to do with the application of the method to ground-based TeV gamma-ray observations. The 
poorer energy resolutiorQ of the ground based instruments {\IS.E\IE ~ 15 - 20%) in comparison 
to the Fermi/LAT resolution (generally |A£'|/£' < 10% above 1 GeV 1.10,1 ) would mean that the 
uncertainty in the dispersion of a single photon (of a few s/Te V) could easily hide any anticipated 
dispersion. 

A number of different approaches exist that are s pec ifically designed for tests of time lags be- 
tween event sequences, such as likelihood methods ll ill and modified cross-correlation functions 
applied to the individual photon events [12]. A particularly attractive and simple algorithm was 
the second approach used in the Fermi analysis of GRB 090510. This method was conceived 
to solve the problem of detecting energy-dependent time lags in statistically limited photon lists, 
and the fundamental idea of such a dispersion cancellation algorithrr|| method was independently 
proposed by Scargle et al. (2008) 1 13] and Ellis et al. (2008) [ 14] - the latter derived actually to 
search for QG signatures from neutrino propagation. We also use this technique as the basis for 
our search methodology, but introduce a different test metric that is better suited to the systematic 
uncertainties associated with a VHE photon dataset, in particular the poorer energy resolution. 

In general, if the expected energy-dispersion is small compared to other relevant variability 
timescales of the astrophysical system under study, its exact functional form is of little impor- 
tance, since the dependency can be treated perturbatively and expressed as the first-order terms 
of a Taylor series (cf. equation[T}. The dispersion cancellation algorithm uses this fact and works 
directly on the time- and energy-tagged events to search for a non-zero parameter t (measured 
here in s/TeV) that optimally cancels any spectral dispersion present in the light-curve. The 
lag-correction, 5f,-, on photon / of energy is given by: 

Sti = -TEf (2) 

where a defines the dominant term of the series expansion for the energy dependency of the 
time lag, usually taken to be the linear expansion term a = 1, or the quadratic term a - 2. The 
dispersion cancellation algorithm cycles through a range of possible values for t, looking for 



^The energy resolution is defined as \hE\IE, where l\E is the difi'erence between the true energy and the analysis- 
reconstructed energy of an event. 

'This name was coined by Scargle et al. (2008) in the context of their particular version of the test, but we will adopt 
it here with greater generality 
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the T* that extremises an appropriate metric (or "cost function"), chosen to quantify the presence 
of spectral lags. An advantage of this approach is that it makes no a priori assumptions on the 
nature of the lightcurve apart from the inevitable hypothesis of simultaneity of emission at the 
source. 

A number of different test metrics have been proposed for the purpose of quantifying the 
spectral lag and finding the optimal dispersion cancellation parameter They all use some kind of 
measure of sharpness of the peak in the burst profile as the value to be maximised in the search 
for T* (see examples in 1IT4I1 . jstl and |13]). Here by "sharpness" of a burst we mean a quantity 



proportional to the gradient of the photon density at the time of the maximum in emission. The 
principle behind the maximum sharpness choice is that whilst the emission of high and low 
energy photons at the source is simultaneous (top left panel of figure [TJ, an energy-dependent 
dispersion introduced during the photon propagation will always skew the overall light-curve. In 
the particular example shown in figure [1] this happens by the delayed arrival of the higher energy 
photons (lower left hand plot), thus skewing and broadening the burst profile as a result. The 
maximally sharp burst configuration will be retrieved when the temporal sequence of events is 
again randomised in energy, corresponding to the exact cancellation of the dispersion. Observe 
that this approach will always give a unique solution for each given dispersion model, because in 
the case of under- or over-corrections, r, the asymmetric effect will either stiU be left present or 
be re-introduced in the opposite direction, and the burst will remain broadened in respect to its 
original width. The cost functions used in [13, ^ serve well to minimise the total inter-photon 
spacing within the entire event sequence, thereby maximising the peak of the lightcurve, but the 
poorer energy resolution of the ground-based instruments limits the efficacy of the method when 
applied to VHE observations of sharp bursts. In this paper we examine an alternative test metric 
based on the Kolmogorov distance between two probability distributions, which better exploits 
the fact that, whilst the energy resolution of an individual photon is far from ideal, the overall 
energy bias of a sample of them is actually 2, AE, ^ 0. 

2.1. The Minimum Kolmogorov Distance Metric 

For a data-set with a sufficient number of photons (a few tens), the event list can be sepa- 
rated into low- and high-energy bands, forming two independent datasets. In the absence of any 
spectral dispersion, the basic assumption that the temporal sequence of events is randomised in 
energy should hold and the profiles (apart from statistical fluctuations plus some arbitrary inten- 
sity scaling that can be eliminated by normalisation) should superpose. If, however, a systematic 
spectral dispersion is present, the profiles of the light-curve will look skewed relative to each 
other (see Figure [T] lower panels). 

Given two random variables X and 7 in M, a simple measure of the difference between their 
respective probability distributions is the Kolmogorov distance Dk, defined as the maximum 



vertical distance between the two cumulative distribution functions (CDFs) II15I1 : 



Dk = sup \Fx{x) - Fy{x)\ (3) 

where Fxix) - prob(X < x) and Fy{x) = prob(Y < x). 

The situation is illustrated in the right-hand plots of figure [T] Assuming that the events are 
generated simultaneously and co-spatially in the source, any energy-dependent dispersion intro- 
duced between the two will show up as an increased D^ between the cumulative distribution 
functions of their arrival times. Therefore, minimising this value will amount to cancelling any 
dispersion present (simultaneously minimising the sharpness of the profile). It is well known 
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Figure 1 : Cartoon of the effect of the energy dependent dispersion on the shape of the low (L) and high (H) energy 
profiles. Observe that the systematic shift on the high-energy curve relative to the low-energy one induces a skew to the 
burst. The panels to the right show the coiTesponding discrepancy in the cumulative distribution function. The maximum 
vertical distance is indicated, con'esponding to the Kolmogorov measure Djc. Observe that Dk tends to fall always in the 
middle of the distribution, near the peak position of the profiles. 



5 



from the properties of the Kolmogorov-Smimov test that the Kolmogorov distance is insensitive 
to the tails of the distributions, where the CDFs converge to the values of and 1, and which 
describe the probability of extreme events lll6ll . In fact, Dk will tend to fall around the central 
regions of the CDF, near to the peaks of the profiles where their accumulated discrepancy is max- 
imum. This is a useful property because it means that the measure naturally attributes a greater 
weight to the most transient parts of the light-curve, whilst being relatively insensitive to outliers. 



3. Performance of the Method 

We now analyse the performance of the method by discussing the four main factors that are 
expected to affect the sensitivity to detect energy-dependent dispersion: burst width (section l3T2l i. 
energy resolution (section l33T l, burst intensity and asymmetry (section [j4t . Before we move on 
to discussing these specific topics, we list the steps for application of the algorithm: 

• select a burst or transient event from the light-curve; 

• split the burst photons into low- and high-energy datasets, this will be a trade-ofi^ such that 
both groups have the largest possible number of events in them, but also that the difference 
between their average energy is as large as possible; 



• build the CDFs for the two distributions (see Appendix A for some discussion on alterna- 
tives ways of representing the lightcurve); 

• adopt a model for the time delay, e.g. linear or quadratic in energy; 

• apply correction to the time-stamp of photons according to equation|2l 

• for each value r of the correction, calculate Dk', 

• the optimum t* is the one which minimises Dk for the range of r tested; 

• assess the uncertainty by simulation of the burst or bootstrap of the data. 

For illustrative purposes, we will only consider in this section the ideal case of an isolated 
Gaussian burst. The superposition of multiple bursts or burst shapes different from Gaussian will 
be discussed when the method is applied to real flare data from the AGN PKS 2155-304 in the 
next section, but do not change the conclusions presented here. For our studies individual burst 
data were simulated using the generalised Gaussian shape fvf^, which can also provide a good 
match to the pulse profiles generally observed from AGNs and GRBs: 



/(f) = Imax exp 



(4) 



where t is the time into the flare, f^av is the time of maximum flux /,„,„, cr,. and cr^ are the 
signal rise (for t < t„,ax) and decay (for t > t„ax) time constants respectively, the "sharpness" 
(peakedness, or kurtosis) of the profile is given by the parameter a- > 0. A low value of k means 
a sharply peaked pulse, a high value a more rounded one, and k - 2 corresponds to a Gaussian 
pulse shape. The rise (fr) and decay itd) times from half to maximum amplitude are found from 
the rise and decay constants using the equation f,-,/ = [ln{2)^^'']crr4. The spectra are assumed to 
take a power law shape of the form dN/dE = kE^, where F is the spectral index and k a flux 
constant. 
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3.1. Energy cuts 

The first step necessary in constructing the CDFs for the analysis means we must decide where 
to place the low- and high-energy boundaries. This choice is made such that the difference in 
the mean energy between the two CDFs can be maximised, while keeping good photon statistics 
in both energy bins for the analysis. We have verified that due to the (usually) steeply-falling 
spectral index of the photon distributions, the analysis is less sensitive to the choice of the low- 
energy boundary, provided this is set comfortably above the threshold energy of the instrument. 
We set here the low energy band to be 0.2 < S < 0.4 TeV. We then searched for an optimal 
high-energy cut window, which will be the more statistically-starved component. Simulations 
were for a Gaussian light curve shape of 120s rise/fall time and a maximum count rate of 3 Hz. 
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Figure 2: Effect of tfie choice of the minimum energy cut for the high energy band (H„,i„) on the accuracy of the 
determined dispersion measure. Simulations are for power law spectra with indices of -2.5 (crosses), -3.0 (squares) and 
-3.5 (triangles). 



Figure |2] shows the results of our analysis on the effect of the choice of the minimum value 
for the high energy cut H^n on the root-mean-square (RMSJI of the reconstructed dispersion 
parameter (t*) for three different values of the spectral index (for F - -2.5, -3.0, -3.5). The 
curves for each spectrum show a slight improvement of the RMS with increasing energy cut 
which quickly plateaus. This is because at too low a minimum energy for the high energy sample 



The root mean square of a distribution X with N events Xj is given by RMS = yS^/W- 
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the difference in the average photon energy of the events gets too low and the high-energy CDF 
becomes indistinguishable from the low energy CDF as they are both dominated by events that 
are not dispersed much. The error in any reconstructed lag correspondingly increases as natural 
statistical fluctuations become the dominant source of uncertainty. Much more noticeable for this 
though, is the steep rise in the RMS for Hmin > 1 TeV for the softest energy spectrum (F = -3.5). 
This occurs when the number of events in the high energy sample drops below ~ 10 - the CDF 
becomes ill-defined with respect to the low energy sample CDF and so no reliable dispersion 
measure can be found. This result gives an idea of the minimum number of events necessary in 
an energy band for the method to be able to work. From this we see that the results of the test 
will be fairly independent of the actual minimum high energy cut, provided the high energy band 
has > 10 events and is at least a factor of two higher in energy than the lower energy sample. We 
take Hmin > 1 TeV for the remainder of the paper, unless specified otherwise. 



3.2. Sensitivity to burst width 

We quantify the sensitivity to the burst width by the term "sensitivity factor", rj, following the 
definition in j^. This quantity is written as the ratio of the expected dispersion magnitude (6t) to 
the width of the transient feature (At): 



This ratio is the main parameter which will quantify the size of the lag that can be probed by the 
method, for a given burst width. 

For this analysis, we simulated 10,000 Gaussian burst profiles of 500 events each, with a low- 
energy threshold of 200 GeV and a spectral index F - -2.5. A dispersion was then introduced 
that varies from 5-200% of the burst width, i.e. from the dispersion being entirely contained 
within the burst to the burst being smeared in time over a period greater than its duration. The 
results are shown in figure |3] where the points correspond to the mean reconstructed r* and the 
error bars are the RMS of that distribution. We see, as expected, that the narrower the width of 
the burst with respect to the introduced delay, the better the delay can be determined. The error 
bars in the plot indicate the 68% confidence intervals (CI) of the reconstructed lag distribution, 
showing that the method can reconstruct the correct value of t and exclude the null hypothesis 
of zero lag at the 99% level up to a value of t; ^ 0.2. This corresponds to a sensitivity limit of a 
lag equal to 20% of the burst width. 



3.3. Sensitivity to energy resolution 

We also included in our sensitivity analysis the effect of the energy resolution (|A£|/£), an im- 
portant consideration in ground-based gamma-ray measurements. This uncertainty will directly 
affect the dispersion correction for individual photons and will thus limit the sensitivity of the 
method. The energy resolution is modelled as 0% (an ideal detector case where the reconstructed 
energy is always the true energy), 10% and 20%, shown as the sub-sets of data in figure[3] There 
is a small systematic trend for the reconstructed lag to be under-estimated as the energy resolu- 
tion gets poorer, but this is very small in comparison to the overall error in the reconstructed r*. 
The under-estimation is expected: the power law nature of the spectrum means that any width 
to the energy resolution will systematically spill more photons to higher reconstructed energies 
than photons to lower reconstructed energies by sheer weight of numbers. This is a well known 
problem in spectral reconstruction of VHE sources. It is possible, with appropriate Monte Carlo 
modeUing or bootstrapping, to compensate for this systematic trend if necessary. 
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Figure 3: Sensitivity of the algorithm to the ratio lag/burst width r] for 0% (open circle), 10% (open square), and 20% 
(open triangle) energy resolution. The results are from sets of 10,000 Monte Carlo simulations of Gaussian profiles, 
containing 500 events each, for a low-energy threshold of 0.2 TeV and spectral index r = -2.5. The low- and high- 
energy bins were defined such that the average energy difference between the two is ~ 1 TeV. 
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3.4. Sensitivity to burst intensity and asymmetry 

The burst intensity is another factor that will affect the sensitivity of the algorithm, since 
it will limit the photon statistics available to construct the CDFs. This is shown in figure |4] 
and was tested over a similarly-generated set of Gaussian profiles as before. Here bursts were 
generated with differing numbers of events, between 50-3000, accounting for different count 
rates (corresponding to /„„v between 1-lOHz) and for 3 different burst widths, with rise/decay 
times varying between 10- 120s. 

For a given burst width, the effect of increasing the number of events in the light curve is 
to reduce the RMS of the recovered dispersion parameter From a certain number of events 
onwards, and depending on the width of the burst, the distribution tends toward a plateau and 
little improvement in the RMS is obtained by further increasing the number of events. As noted 
earlier, the sharper the burst, the earlier this plateau is reached. Finally, we have also tested for 
any effects due to profile asymmetry by maintaining the total burst width and varying the ratio of 
rise/decay time of the flare. The results plotted in Figure |4]show that the method is not affected 
by any intrinsic burst asymmetry, but only by its overall width. 
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Figure 4: Sensitivity of tlie metliod in relation to the width of the burst and the number of events in it. The labels in the 
key define respectively the rise and fall times of the profile. The results are from MC simulations of 10,000 bursts with a 
maximum count rate between 1-lOHz. 



3.5. Sensitivity to energy spectrum 

The observation of a single burst or flare is not going to provide definitive evidence (or refu- 
tation) of energy dependent dispersion due to QG. Instead, a number of sources demonstrating 
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a consistent behaviour for a range of redshifts will be necessary to be able to confidently deter- 
mine if such an important effect exists, and to disentangle it from source-intrinsic lags. Even 
if the intrinsic spectrum for a given source type is identical between objects, the interaction of 
the gamma-rays on the diffuse extra-galactic background radiation will lead to a softening of 
the observed spectral index with increasing redshift. The number of very high energy events is 
intimately related to the energy spectrum; to quantify the systematic uncertainty introduced by 
this effect in the estimate of t*, we simulated 10,000 Gaussian shaped lightcurves, with 120 s 
rise and fall times and a maximum count rate of 3Hz, for a range of indices T between -2.5 and 
-3.5. The RMS of the recovered r* is plotted in figure |5] demonstrating an approximately linear 
deteoriation on the determination of r* as the spectrum softens. This effect is also easy to un- 
derstand as the softening of the spectrum represents a depletion of high-energy photons from the 
high-energy CDF. 



24 
22 
20 
18 

'in 

12 
10 

8 

6 

2.4 2.6 2.8 3 3.2 3.4 3.6 

Spectral Index 

Figure 5: The effect of the spectral index on the uncertainty in the recovered dispersion. The simulated bursts have a 
width of 120 s, meaning that the recovered dispersion's RMS varies from 0.5-20% of the burst width, for an index T 
going from -2.5 to -3.5. 




3.6. Burst analysis window 

Until now we have treated simulated isolated bursts, for which we are confident it is straight- 
forward to find a reasonable signal-to-noise ratio above background from which we can define 
the burst to start/end. When analysing transient events within a real hght-curve, as will be done 
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in the next section, it is important to consider the effects of confusion and under-sampling of in- 
dividual bursts. If the burst is adjacent to other structures within a complex lightcurve it might be 
difficult to define with precision its start and end times, and a superposition of different features 
might be unavoidable. In particular, the highest energy, most-lagged events could then fall out- 
side an inappropriately chosen analysis window, thus affecting the profile reconstruction. Also, 
if the burst is at the edge of an observation run, data could be missing for part of the flare, this 
loss of information will also be energy-dependent if there are lags in the light-curve and this is 
likely to affect the performance of the reconstruction. 

To test for these effects and assess if a proper reconstruction of the original profiles of lagged 
light-curves is still possible within our framework, we performed two sets of MC simulations, 
for which we generated two groups of 10,000 Gaussian bursts with 500 events each, a spectral 
index of F = -2.5 and an energy resolution of 20% was used. 

For the first set, represented in figure |6] the analysis considered a series of windows around 
the peak position of the burst of widths equal to 1, 2, 3 and 5 times the combined rise/fall time 
(tr,d) of the burst (where these relate to the time for full-width-half-maximum t,- + td - tfWHM), to 
simulate different degrees of under-sampling. In this case, a so-called "transparent window" was 
used. This means that the CDFs are built only with the events that fall within the time window 
boundaries after the dispersion cancellation has been applied, but for each different value of t 
events are allowed to pass into and out of the window's boundaries - updating the CDF at each 
new step of the algorithm. 

The result is that a very narrow window around the burst affects the accuracy of the recon- 
struction, increasing the RMS by up to 20% when only the FWHM around the burst peak is used 
to build the CDF. This degrading effect can be understood as a consequence of an ill-defined 
shape for the CDF. The effect is present for all the range of sensitivity factors tested, being more 
pronounced for smaller This suggests that in using the method one should attempt to include 
as much of the burst as possible into the analysis, in order to include the most possible infor- 
mation on the profile shape for the CDF comparisons. An arbitrary choice of a narrow window 
about the peak of the burst to artificially reduce rj does not improve the results, due to a loss of 
information in the CDFs on the shape of the wings of the distribution. Of course this observation 
is no prescription for the analysis. Ideally the time window around the burst should be at least 
^tfWHM, but in the case that the burst is confused with other features the analyser should simply 
be aware of this degrading factor when determining the confidence intervals for r. 

This also demonstrates that whilst individual sub-flare features may be resolved, they can still 
potentially influence each other If a train of bursts is too closely aligned with respect to each 
other for an individual analysis to be conducted, then our simulations show that the critical factor 
is to account for the rise/fall parts of the profile of the first/last burst respectively. In such cases, 
the best approach to the analysis is to consider the train of bursts together, rather than trying to 
split the bursts into ill-sampled individual features. 

As a final note, in the same way that events pertaining to the burst can be selected out of the 
analysis window, events not pertaining to the burst can also contaminate the analysis during the 
cancellation procedure, for instance from background events or from superposing bursts closely 
aligned in time. This will produce similar kinds of effects as the case treated in figure |6] As 
before, a compromise has to be found (if necessary via dedicated simulations with a configuration 
as similar as possible to that of the real lightcurve) between reducing the contamination and 
sampling well the profile of each individual burst. 

There is a second class of situations when a burst will be under-sampled at the detection level 
rather than in the analysis procedure, for instance when the burst occurs close to the start/end 
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Figure 6: Sensitivity of the Kolmogorov distance metliod in relation to the size of a "transparent" window used to 
reconstruct the CDFs from the burst profile. The labels in the key define different data sets with different sensitivity 
factors = 0.5 (star), 0.2 (triangle), 0.1 (circle). The results are from MC simulations of 10,000 bursts generated from a 
Gaussian shape and an associated energy resolution for each event of ~ 20%. 



13 



times of an observation run. In this case, events are lost permanently. To simulate this effect we 
introduce a so-called "opaque window" in the analysis, where events that are initially out of the 
analysis window for a r = will not be admitted in as t changes during the steps of the cancella- 
tion algorithm. The results are shown in Figure Q In this case, the fact that we permanently lose 
a high proportion of high energy events (because they are more spread or dispersed then the low 
energy ones) means that not only will the RMS be worsened, we preferentially lose the events 
that will most accurately recover the correct dispersion. 

The three different data groups represented in Figure Q are for sensitivity factors 77 equal to 
0.5, 0.2 and 0.1. Note that the case of longer duration bursts is the most affected, simply because 
in this case more high-energy events are permanently lost from the burst window, relative to the 
low-77 case. Within each dataset, points 1-5 indicate the size of the window in units of tfWHM- 
Here, the permanent loss of information about the most-lagged events mean that the true lag 
T* is reconstructed wrongly the closer the burst maximum is to the start/end of the observation 
window. Therefore, as before, the conclusion is that windows as wide as ~ 3tfWHM around the 
flare peak give the best compromise between burst width and information content. Another con- 
clusion from this analysis is that the search for lags in a flare for which a large portion is missing 
from the burst (more than 1/4 of the the total number of events) is certainly not recommended, 
so observation strategies that maximise the amount of on-source time in lightcurve monitoring 
(such as fis"]) are definitely encouraged. 

In any case, as before, a complete simulation, or bootstrap, of the observed light-curve is the 
best way to assess the correct RMS from the combined effect of all cases discussed in this section. 
An example of this approach to the analysis of real data-sets is given next. 

4. Application to PKS 2155-304 flare simulations 

In the previous section we have discussed all the principal factors contributing to the uncer- 
tainty in the reconstructed lag. We have illustrated those by the simple case of isolated Gaussian 
bursts. We are now in a position to test the efficacy of the method to recover a dispersion in 
some realistic lightcurves. To do so, we move from the simple tests on Gaussian profiles to work 
on simulated datasets based on fitted profiles for the large flare of PKS 2155-304 from 2006 
101. For consistency with previous work done on the PKS 2155 data set, and to enable people 
to reproduce our work, we use exactly the same profile fits presented in the original H.E.S.S. 
publication, instead of searching for and separating the individual bursts ourselves (step 1 of the 
analysis procedure). This introduces a binned and parameterised aspect to the method, for that is 
how the lightcurve features were identified in [7]. To keep the method truly non-parametric and 
unbinned, an approach such as the Bayesian Blocks Il9ll algorithm to search for the time window 
cuts could be applied and is recommended, but it does not change the results we present here as 
this exceptional flare of PKS2155-304 is well resolved. 

There are 5 prominent flaring events (labelled BFl-5) noted for this lightcurve, reproduced 
in figure [8] and the relevant parameters for the generalised Gaussian fits are reproduced in ta- 
ble [1] The simulated event times are generated by random draws from a distribution described 
by Equation m each flare summed to give the total lightcurve. To each event time, an energy 
value is then randomly attributed from a power law distribution, with Ey > 120 GeV. As there 
was no evidence for any spectral variability (at the AF > 0.2 level |i7J) during the night, for 
simplicity we adopted a simple power law distribution with a spectral index of F = -3.5 in our 
simulations; changing the model to the broken power-law fit given in fT] makes no difference to 
the general conclusions discussed here. The error in the energy reconstruction of a single event 

14 



CO 
CO 

-*—> 

CC 
T3 




0.5 1 1.5 

reconstructed lag [arb] 

ri = 0.1 0.2 



2.5 



0.5 



Figure 7: Sensitivity of the Kolmogorov distance method in relation to the size of the "opaque" window used to construct 
the CDFs from the burst profile. The different datasets correspond to = 0. 1 (circles), 0.2 (triangles), 0.5 (stars) respec- 
tively. The number 1-5 within each dataset are for windows of 1-5 (tfwHM) widths, respectively. The results are from 
MC simulations of 10,000 bursts generated from a Gaussian shape, and an associated energy resolution to each event of 
~ 20%. The actual value of the induced lag is indicated by the dotted vertical line. 
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is dominated by systematic uncertainties and is estimated to be of the order of 15% throughout 
the entire energy range. In reality the energy resolution is a function of energy, improving for 
higher energies, and so this value can be taken as a worst case scenario. Then, to simulate the 
energy dependent dispersion, a systematic delay t was applied to each photon's true energy and 
the recovery procedure was carried out based on the instrumentally smeared energy values. 
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Table 1: Parameters used for the generalised Gaussian fit to the PKS 2155-304 flare simulations, based on the original 
H.E.S.S. analysis results |7]. The third column (Max Rate) refers to the maximum count rate of each burst, corresponding 
to its peak flux at time t,„a^. The parameters cr,- and crj are the rise and decay time constants of each burst and k a measure 
of the sharpness of the peak (see text for details). 

In a real-case analysis like the one shown here, there are two non-trivial steps (numbers 1 and 
2 of the list shown in section 3) that must be considered: (i) the choice of the analysis window 
around each burst and (ii) the choice of the energy boundaries to construct the CDFs. Figure |9] 
shows the results of our analysis on the effect of the choice of the high energy cut on the RMS 
of the reconstructed dispersion parameter for the individual flares. The improvement in the RMS 
as the cut moves away from the soft energy band is notable. It is followed by the presence of 
an optimal plateau around and above 1 TeV and worsening RMS above 2 TeV due to a loss of 
event statistics. All this reproduces what was seen in the ideal flare shape case of section lTTI and 
shows the choice of 1 TeV for H,ni„ to be a good one. The uncertainty in the reconstructed lags 
(in s/TeV) were determined from Monte Carlo simulations performed for each individual burst. 
The top panel of figure [10] shows an example of values of the Kolmogorov distance Dk for each 
different dispersion parameter r* tested in the analysis. The middle panel shows the distribution 
of Dk versus t* for 10,000 realisations of BF2, from which confidence intervals for the lag were 
derived, as shown in the lower panel histogram. 

4.1. Testing the recovery of a known induced dispersion 

To demonstrate the efficacy for recovering a dispersion, for each simulated lightcurve in fig- 
ure [TT] we introduced an artificial dispersion between -100 < t < 100 s/TeV, which we aimed to 
recover with the dispersion cancellation algorithm and minimisation of the Kolmogorov metric. 
The algorithm was applied to each of the five major burst features in the dataset, BF 1-5, gen- 
erating five sets of independent measurements. Whilst the RMS of the recovered dispersion r* 
leaves uncertainty in the true dispersion t, the mean of the recovered dispersion is an accurate 
reflection of the true dispersion. This shows that the measurement of multiple flares from a single 
object could be used in combination to give a more accurate estimate of any induced dispersion, 
since they should be the same for all flares if they have an origin in propagation effects. Also, the 
accuracy of the recovery over a large range of the parameter space demonstrates that a number 
of objects at different redshifts could then accurately trace a distance dependent (i.e. propagation 
induced) dispersion. 
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Time [s] 

Figure 8: The pai'ent population light curve of the MJD 53944 PKS 2155-304 big flare event simulation. The individual 
bursts BFl-5 (see tabl^are marked by the thin lines and the dashed line shows the constant signal (~ 10% of the 
maximum count rate) onto which the bursts are superimposed. Values are renormalised from the flux values in [7] to 
count rates here. The heavy line denotes the cumulative light curve. The grey shaded regions mark out the location and 
extent of the l-ti-j windows for the bursts, the black bands correspond to the data gaps due to observing run transitions 
(see section [32). 



17 




1 .2 1 .4 1 .6 
energy [TeV] 



2.2 



Figure 9: Effect of the choice of the energy cut for the high energy band on the accuracy of the determined dispersion 
measiue based on Monte Carlo simulations of the burst profiles BFl-5 (see text for details). 
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Figure 10: Top panel: Kolmogorov distance profiles for diff'erent values of r from the the search of energy-dependent 
dispersion in the bursts BFl-5 of PKS 2155-304. The analyses were perfonned in steps of 1 s/TeV. Middle panel: 
2-D histogram of the Kolmogorov distance calculated for each step of r in 10000 simulated lightcurves for the profile 
of BF2, with no dispersion introduced. A minimum is always reached with some variance around s/TeV. Bottom 
panel: Histogram of r* for each of the minimum values found in the simulations of BF2 shown in the middle panel, 
centred on 0. The standard deviation of a Gaussian fit to this distribution is used to estimate the limiting accuracy for 
reconstruction of r* . 
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Figure 1 1 : The accuracy to which the method can recover a fixed dispersion parameter introduced into the light-curves 
representing each of the sub-flares of PKS 2155-304. Each point is the average dispersion estimated from 10,000 sim- 
ulated lightcurves and the errors represent the RMS of the recovered dispersion parameter, r*. Values corresponding to 
the dift'erent flares are slightly offset on the x-axis for visual clarity. The dashed line is a guide to the eye of y = x, not a 
fit to the data. The dot-dashed line represents the expected linear dispersion for a Planck mass scale of quantum gravity 
from a source at the redshift of PKS 2155-304. 
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4.2. Estimating the limits to QG measurements from the Kolmogorov metric 

It is known from previous studies of this dataset that there is no spectral dispersion present in 
the PKS 2155-304 big flare night lightcurve that is large enough to be detectable within current 
instrument's sensitivities. We therefore limit ourselves to using the simulations of the dataset 
as a test-bed for assessing the sensitivity of the method in a real-life case. The Eqg sensitivity 
limits that could be placed from the application of the algorithm to the big flares from PKS 2155- 
304 are presented in table |2] We adopted a linear relation between the lag and energy of the 
photon (a = 1 in equation |2l) for the Eqc limit estimates in the final column. The mean energy 
difference {AE} between the low- and high-energy profiles is also given. Note that the sensitivity 
limits obtained here are comparable with the most constraining results achieved to date from 
more complex analyses of this particular flare [2]. This result shows that this simple method has 
the potential to probe QG effects to at least the same levels achieved by the current tests using 
AGN, with the advantage that, by separating the light-curve into multiple bursts, from a single 
dataset we can derive multiple independent measures of the same quantity. 
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1.42 


+24.1 


> 0.15 
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1.23 


+ 16.4 


> 0.19 
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1.40 


±18.6 


> 0.2 


BF4 


1.25 


±24.4 


> 0.13 


BF5 


1.29 


±24.5 


> 0.14 



Table 2: Results for each of the PKS2155 big flare night sub-flares (AE) is the mean energy difference between the low- 
and high-energy CDF. o",n„ coiTesponds to the standard deviation of a Gaussian fit to the histogram of r* from 10,000 
simulated lightcurves. The sensitivity limits on Eqq correspond to the 95% confidence levels from the uncertainty in 
the recovered r*. Note that these values are not real estimates of Eqq limits from the PKS 2155-304 observations, but 
sensitivity estimates to searches in a simulated dataset. 

It is apparent in figure [TT]that any dispersion |t| < 50s/TeV is unlikely to be confidently de- 
termined from an individual flare observed with the current generation of Cherenkov telescopes. 
The exceptional flares described here each have too small a value of 77 to resolve such a signature 
individually, the even longer typical widths of observed VHE AGN flares only serves to reinforce 
this point. Nevertheless, the fact that lags can be accurately reconstructed over a large range of 
T, demonstrated by the correct mean value of the reconstructed lags, supports the idea that the 
combined results of a number of flares from an individual source could improve the Eqg limits 
considerably for each given object. 

5. Discussion and Conclusions 

We presented here an unbinned, non-parametric method of testing for energy dependent dis- 
persion in a lightcurve that is sufliciently robust to work under the constraints of scarce counting 
statistics and the modest energy resolution expected for VHE gamma-ray data analysis, the latter 
issue having been noted as a limiting factor in other simple dispersion cancellation methods, eg 
fisll . The applicability of the method is based on the fact that a random distribution of events 
in energy will give rise to indistinguishable time profiles (or CDFs) that can thus be directly 
compared using some kind of statistical metric, among which the Kolmogorov metric was found 

21 



by us to be best suited to our purposes. When used in conjunction with unbinned algorithms 
for identifying variability features like flares in lightcurves, such as for example the Bayesian 
Blocks algorithm lll9ll . the analysis chain would be an entirely non-parametric way to search for 
energy dependent dispersion (though such an analysis is beyond the scope of this rather specific 
discussion). 

We have discussed in detail the factors contributing to the uncertainty in the determination 
of a lag, and the method was subsequently applied to the challenging case study of looking for 
a small magnitude dispersion expected from a specific QG-induced LIV, following the form of 
equation [1] It is akeady known that there is no dispersion present in the PKS 2155-304 big flare 
night lightcurve that is any larger than the limits achievable by our method, presented in table |2] 
yl, therefore we did not attempt to derive further limits using that exemplary dataset. It should 
be noted, however, that the limits for each of the simulated sub-flares are at least a factor of 
two better than the cross-correlation method of fT] and comparable to the much more complex 
analysis of |i2J; but while those two methods used the entire night's lightcurve as a single dataset 
to derive their limits this method has the advantage of being able to treat distinct bursts within 
the lightcurve as independent tests. Analysing a larger number of individual sub-flare features 
with the method presented here, and posteriorly combining a statistical sample of flares to obtain 
a single estimate (or limit) on dispersion could correspondingly make for yet more constraining 
limits on LIV. 

In section 13.51 we found the intuitive result that the uncertainty in any measured dispersion 
scales inversely with the hardness of the energy spectrum. This also leads to the slirfitly counter- 
intuitive consequence that a nearer AGN such as Mrk421 in a high state (eg 

m, M) could 

provide as constraining a limit, if not more so, than a more distant object (like PKS 2155-304) if it 
showed similar variability timescales. This is solely due to the harder observed spectra, meaning 
an increased number of high energy photons being detectable from the nearer object (due to 
less absorption by the extragalactic background light). This could give a smaller uncertainty 
on the dispersion even if the expected magnitude of the delays are smaller At a redshift of 
z ~ 0.03 for Mrk421, the expected dispersion would be t ~ 1 s/TeV for Planck mass scale 
quantum gravity; with a next generation observatory such as CTA [22] the number of photons 
above lOTeV in minute scale flares could potentially push the expected sensitivity t] an order 
of magnitude beyond what is achievable today with AGN studies. This means that such AGN 
studies could perhaps surpass even the current Fermi limits on GRB 090510 [6] making for a 
useful independent confirmation of that result. The combined information from different source 
classes and over a sufficient range in redshift would then enable to distinguish between intrinsic 
source and external propagation induced dispersion effects. 

The performance on the method has been tested on a specific case, but this method is, of 
course, not just limited to searches for a linear expansion term and VHE gamma-ray observations 
provide the highest energy photon datasets for probing the quadratic term of equation[T](e.g., ||2l). 
The simple isotropic dispersion scheme of equation[T|also neglects a number of aspects of poten- 
tial Lorentz violation eff'ects, such as birefringence and direction dependent vacuum co-efficients 
for Lorentz invariance (see eg equation 145 of [23]) which could make up interesting follow-up 
studies with either this or next generation.gamma-ray telescope datasets. Neither is our method 
just limited to ground based VHE gamma-ray datasets, the extended baseline for the Fermi GRB 
datasets can more than compensate for the lower photon energies (with regard to the linear term, 
l^b), and tests for birefringency effects (which instead cause a spread in the wave packets 1231') 
have been made using INTEGRAL observations of the highly polarised GRB 0412 19A. Unfor- 
tunately it is not currently possible to determine the polarisation of gamma-rays from ground 
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based obervations, so this case has not be deah with here, but broadening of the Hghtcurve can 
also arise in certain models that predict a stochastic, super-luminal aspect to the refractive index 
(eg II24II ). though such tests may be better aided with a metric that is optimal in testing the change 
in kurtosis of the lightcurve, rather than the skewness, as used here. 

Most importantly of all, the method for testing for the presence of a dispersion is not limited 
to just testing forLorentz invariance violation effects. Indeed any physical model that introduces 
a width or a delay in the emission, such as acceleration within the emission region (eg (t^]) or 
due to cascading on the intergalactic magnetic field (eg 1,261 ). can be adopted as the model for 
the time delay correction. These are all tests to be pursued in future work and not in the scope of 
this paper, which is to present the method, its benefits and limitations. 
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Appendix A. Light Curve Representation 

In order to apply the algorithm, it is necessary to define how to construct the CDFs from the 
low and high-energy event sequences of the burst under study (step 3 of the analysis procedure). 
Given that the Kolmogorov distance is a metric for probability distributions, the event sequence 
must be normalised. Since the dataset is composed of time- and energy-tagged events, the can- 
cellation will be applied to every photon individually so that none of the available information 
is left out of the analysis. The simplest choice for representing the data is therefore to construct 
empirical CDFs for both the low- and high-energy profiles as step functions from the original 
event sequence, according to the following rule: 

CDF : F{ti) = i/N, (A.l) 

where f, is the time of the i'* event in the sequence, and is the total number of events in the 
sequence. In this construction, the height of each step is constant and equal to ' (the CDF is 
defined between and 1), and the length of each step equals the waiting time between events in 
the sequence. All the timing information of the temporal sequence is thus explicitly preserved in 
this representation. 

A different representation of the dataset was proposed in Scargle et al. (2008) [13], and can be 
used as an alternative way of constructing the CDF. In this representation, the dataset is tesselated 
so that the photon sequence is represented by a series of cells of width dt, constructed around 
each event /. A cell density is then defined by the rule x, = 1 /dtt, which can be interpreted as the 
instantaneous rate of the process at time f,, which is later normalised into a discrete probability 
distribution: pi = jc,/ 2 Xi. The CDF in this case would be: 

CDF : F{ti) = 2 Pi, (A.2) 

For the application of the Kolmogorov distance metric, we found that the first representation 
in equation lA.il is more appropriate. This is because the magnitude of the cell densities repre- 
sentation can be dominated by spikes resulting from very small inter-event times in some cells, 
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Figure A. 12: Choice of the light-curve representation. The panels on the left show the cell density representation for the 
low- and high-energy components of flare BFl of PKS 2155-304 (see|4]for nomenclature of the flares). The right panels 
show the correspondent CDFs for the two Ught-curve representations discussed. Note that the raw events representation 
shows considerably less "raggedness" in the CDFs than the cell density one, and is therefore more appropriate for 
calculating the Kolmogorov distance cost-function. 
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originating from the noise of the Poisson process, which will introduce excessive "raggedness" 
in the CDF representation. This can be seen in the right panel of figure lAT2l where we compare 
the low- and high-energy CDFs from a simulated burst profile based on the burst BFl of the VHE 
flare of PKS 2155-304 observed with H.E.S.S. In this case, both profiles superpose, as there is no 
significant dispersion present, but the cell density representation results in additional fluctuations 
in the constructed CDFs. A way to circumvent this problem within the cell representation is to 
adopt a logarithmic scale for the density - for example x, - log( 1 / dti) - which better recovers 
the shape of the profile. 
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